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Abstract 

Using a recently developed method for the inclusion of fluctuation in the 
BUU dynamics, we study the self-consistent propagation of inherent thermal 
noise of unstable nuclear matter. The large time behaviour of the evolving 
system exhibits synergism between fluctuation and non-linearities in a uni- 
versal manner which manifest in the appearance of macroscopic structure in 

the average description. 
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In recent years a considerable interest has grown in the understanding of the role of fluc- 
tuation in the dynamics of multifragmentation. In the decompression phase of the collision, 
the density of the composite nuclear system eventually attains such a low value so that it 
crosses the boundary of the so called spinodal zone where a small perturbation can grow 
rapidly until the fragments form. The fluctuation certainly plays a significant role in the 
process of fragmentation. To unravel some basic features of such delicate process simple 
model study seems encouraging. For this purpose, different theoretical approaches based on 
Landau [|],^ or Boltzmann-Uhling-Ulhenbeck (BUU) P] equation have been utilized. 

Let us concentrate on the BUU approach to the problem. In this model, the description 
of the process is governed by the evolution of reduced one body density distribution. In the 
mean field approximation, the fluctuations arising from the stochastic part of the collision 
term are ignored and evolution of a single average effective density is considered only. In this 
description the bundle of trajectories are confined in a certain region of the associated phase 
space. The presence of collision term in the equation corresponds to the dissipative effect 
in the evolution. When the system reaches the point of instability where any fluctuation, 
however, small it is, amplifies indefinitely so that trajectories may bifurcate and drive the 
system to some new fixed point. In other words, different manifestation of the system are 
now become accessible dynamically. Therefore, the mean trajectory description has been 
lost. To encounter such a situation, the general strategy is to study the evolution of a 
collection of initially prepared samples. Due to the very presence of the stochasticity in 
the process their fates are different. To make this picture consistent, a usual method is to 
introduce a fluctuation term in the equation of motion, so that the nature of the solution 
become probabilistic. Within the framework of BUU equation this prescription is properly 
incorporated and the resultant evolution equation for phase space density is commonly 
known as the Boltzmann-Langevin (BL) equation 

Guided by the same physical principle, several computational schemes have been pro- 
posed to incorporate fluctuation in the BUU equation 1^,^. In this context, we want to 
mention a very recent development (see Ref. where a novel method of simulation is pro- 
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posed to incorporate fluctuations in the BUU equation. The occupancy factor n(r, p) at 
a particular elementary cell As of size ArAp//i^ (D the dimension of the physical space) 
around a point (r, p) in the associated phase space of the one body description of the system 
is considered to be discrete in nature so that it takes the value 1 (say success) or (failure) 
otherwise. The basic transition probability (say, 12 1'2') within a small interval of time is 
given by T(12 1'2') = to'i2^i'2'AtAs^nin2ni''n.2' where Ci;i2-^i'2'As'* is the basic transition 
rate. For a particular transition, every sample is tried with a probability of success cuAtAs^. 
However due to the presence of four occupancy factors in the above relation a successful 
transition occur when an appropriate configuration of particle and hole states appeared in a 
sample. By this way enough fluctuation in sample space can be generated and at the same 
time, the relevant statistical criteria i.e. = {n)i{n)i can be fulfilled exactly for every cell 



Another essential aspect of the simulation is to generate the solution of the Vlasov part 
of the BUU equation in the presence of the mean field potential. Due to the discrete nature 
of the phase space occupancy n(r,p), the usual procedure of matrix method |^,^ cannot 
be applied here properly. To consider Vlasov propagation of phase space distribution we 
follow the Nordheim's approach, closely related to the well known test particle method of 
solving Vlasov equation. A filled unit phase cell may be considered as a particle. The time 
evolution of these particles can be described by the standard leap-frog algorithm in the 
following manner. 



where (X^(t), P^(t)) represent the kth component of position and momentum of the ith 
particle at time t. It is important to note that unlike to the test particle method the position 
and also momentum of the particle are represented by the grid points here, therefore, its 
evolution with time take places in discrete steps. However to provide a reasonably good 
energy conserving evolution, we incorporate a modification in the leap-frog routine. The 



1. 



Xl{t + At)=Xl{t) + -PlAt- 

m 

Plit + At) = Plit) + VfcW(p(r, t))At 



(1) 



(2) 
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error comes mainly due to the rounding off the quantities XI/ Xd or P^/pd where Xd and pd are 
the co-ordinate and momentum grid spacing. The accumulated error that arise in momentum 
due to the incoming particles at a particular spatial cell (say I) is given as Ap^ = J2i ^P\- Here 
5p\ is the residual part of PI so that PJ. = pd x nearest integer oi{Pl/pd) + 5p\. Therefore, 
to reduce this error substantially, one may assume that the momentum of all the particles 

pi Ao^ 

in a cell / are modified as pd x nearest integer of(-^) + Here Ni is the total number of 
incoming particles in cell /. To check the reliability of this method we consider the evolution 
of the breathing mode of a spherical nucleus in 2-dimension having small skin of about 2 
fm. We observe two interesting points: 1) the initial spherical symmetry in momentum 
and co-ordinate space are maintained exactly and 2) no two particles appear simultaneously 
in the same phase cell throughout the evolution. However, it is to be noted that when 
fluctuation in phase cell is taken into account by allowing the the collisions among the 
particles the above mentioned features are gradually lost. In this situation, the smearing of 
the distribution becomes necessary. This is incorporated by simply moving a particle in the 
neighbouring sites in momentum space if it tries to cell in the process of evolution 

which is already filled. In this simulation, the values of Xd and pd are taken as 1 fm. and 
40~60 MeV/c respectively. The total energy per particle is found to be conserved up to 
95% even at time t = 100 fm/c. 

Using the above mentioned simulation procedure the present study revisits the problem of 
nuclear matter in spinodal zone. This work mainly concerned about the large time behaviour 
of the system and the non-linearities that appear in the process of evolution. To make this 
study consistent with earlier investigation we adopt the scenario that was studied rigorously 
by several authors namely, a gas of fermions situated inside a two dimensional torus. 
For the effective one body field at two dimensional grid point (x,y) we employ a simplified 
Skyrme interaction U{x,y) = + B{£^f with A = -100.3 MeV, B = A8 MeV, 

the saturation density po = 0.55 fm~^. To evaluate effective potential on the lattice an 
averaging over 2-dimensional Gaussian function having width 0.87 is done. The in-medium 
cross section is taken as 2.4 fm and Fermi-momentum Pf = 260 MeV/c. For simulation 
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we consider the square lattice of size (21x21) in co-ordinate space having width = 1 
fm., (31x31) in momentum space with = 40 MeV/c and the number of samples AT^ for 
parallel runs are taken to be 100. To ensure that our system is initially situated inside the 
spinodal zone we prepare the samples having uniform (in space) average density (p) (or the 
number density (A^))=0.5po (^o) and momentum distribution same as that of a fermi gas 
at temperature T = 3 MeV. 

At a given temperature T ^ there exists thermal fluctuation in the initial state. To 
incorporate this important feature within our simulation procedure we prepare the sam- 
ples with the assumption that the distribution of is given by a Gaussian function with 
mean (A^) where as the sample average value of occupancy {n)i is given by the fermi dis- 
tribution of appropriate temperature T. The variance afj- satisfies the standard relation 
a% = J2NJ2<^i = J2NJ2{^)i{^)i- To maintain spatial uniformity in initial state, the value 
of A^ or the momentum distribution is considered to be the same at every spatial cell of 
a particular sample. Inspite of the fact that at T = all collisions are Pauli blocked, the 
symmetry of the system remains intact even by BL treatment used here or the same used 
in Ref. [3]. It is to be noted that in this case a% = i.e. all samples are identical. The 
presence of the thermal fluctuation in the initial state is in fact, the source of disturbance 
which generates irregularities subsequently in co-ordinate space and may amplify due to 
the action of the mean field. Hence, to make the picture consistent, the preparation of the 
samples in the aforesaid manner is therefore an essential part of our simulation procedure. 

Let us now turn to the investigation of the Boltzmann-Langevin dynamics on lattice. The 
normal mode analysis of the evolving density will allow us to study the interplay between 
instabilities and fluctuations. Following the Ref. we introduce two point correlation 
function to the associated density fluctuation as 

Mt) = {\6pum = f dr' f dre'^-\^'-^\{6pir')5pir)) (3) 

For theoretical analysis of the time evolution of the quantity cr^ the starting point is to 
linearise the Vlasov equation and derive the equation for normal modes (for example see Ref. 
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P). In the presence of fluctuation in the transport process the evolution of the amplitude of 
the Fourier components (or in other words, that of the collective mode u) can be represented 
by the following Langevin equation 

= tu^A^it) + B^it) (4) 

where Ai, is the two component amplitude matrix (with elements ^4+ and A~) of the Fourier 
components of the density fluctuation 6fi, so that 6fi, = A+e*"^"* + A~e~^'^'"^. By{t) is the 
kicking term of the Langevin equation. The associated correlation or diffusion matrix for 
Gaussian Markov process can be given by a Hermition matrix with components 

{m)m')) = 'D^5{t~t') (5) 

where i (j) runs for states '+' or 

Let us now concentrate on the situation of instability. For infinite system one can find a 
range of v , for which all the modes are unstable, the associated the frequencies ujy{= —iVty) 
are purely an imaginary numbers and beyond that range all modes are stable. Therefore, in 
course of time the fate of the system is dictated by the evolution of few such most unstable 
modes. The probabilistic evolution of A^ (see eq.(4)) can be described as well in terms of 
Fokker-Planck equation. The evolution equations for first and second moment related to it 
can be written as DDI 



dt 

da'j{t)i 



{A.m (6) 



-n^^ 

nyj{t) + ni"^ar{t) + 2vi^ (7) 



dt 

where a^^ is a symmetric matrix. The initial condition i.e at t= 0, da^/dt = {a^, = Y.i>j f^u) 
satisfies if T>^^ = Pj7^ = — These equations are the same as those derived in Ref [§, 
the solutions of which provide a good fit to the simulation result |T0[. 

Within the framework of Langevin or FP equation, further study can be directed to 
the non-linear regime of the dynamics by introducing non-linearity in the drift term. Since 
—dUy/dAjj correspond the drift coefficient of FP equation this can be done by adding a 
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quartic term to the associated quadratic potential. This stochastic model of non-linear 
Langevin equation, is rich in physics capable of explaining the general behaviour of the 
system in the presence of instability far from the equilibrium and onset of macroscopic 



structure The presence of the quartic term prevents an indefinite fall of a system and 
drive it towards any one of the two newly appeared stable fixed points which correspond 
to the minima of the dynamically generated cusps. The co-operative behaviour or the 
synergism of nonlinearity and random force simulate this typical feature i.e. slowing down 
the fluctuation near dynamical phase transition. The fluctuation deviate gradually from 
its initial Gaussian profile as the system approaches to the second or scaling regime. The 
modified equation for second moment with positive value of Q is given by 

|crj-+(t) = 2{n, - g^at^{t)}at^ + (8) 

which provides a scaling solution 

-t^{t) = ; r = ^{atHO) + '^)e'^^^ (9) 

with > 0. However, in the derivation of Eq. (8) a physical approximation i.e. A'^'^^{t) = 
(T++(t)yl++(t) is made [O. The time needed to reach the saturation for a unstable mode u 



is given by 



tl = ^log[^W:^0) + ^}r' (10) 



In Fig. 1. the open diamonds represents the time evolution of for unstable modes 
with values of k^^ ranging from 0.3 to 0.6 with ky = 0. For direct comparison to earlier 
results P we plot the function L^a^ (see Fig. 1.) so that the unit becomes fm~^. The 
solution of Eqs. (7) (valid for linear regime of the dynamics) are shown by dotted curves. 
The solid curves represent the solution of non-linear generalization so that the evolution of 
++ component is given by Eq. (9). For these solutions we use the initial condition i.e. at 
t = 0, {A)y = cTj, = da^/dt = 0. The general agreement of our fit to the simulation result 
is quite good. The saturation time for different modes are shown in the figure. The 



most unstable mode i.e. the mode having minimum growth constant tk (= f^^^) saturate 
earhest and the time gradualy increses as the increases. For some cases, the hnear 
regime of the dynamics persists even at time t = 120 fm/c. Let us consider the evolution 
of density p{x,y;t) in co-ordinate space. In earlier time up to ^ 75fm/c when spatial 
average of \f{{S p)^) / {p) reaches the value 1 we observe a considerable fluctuation in p in 
the individual samples. However, the sample average density i.e. {p{x,y;t)) is found to 
be almost uniform with ~ 3% fluctuation from its initial value of 0.5 which is essentially 
the mean (or single) trajectory result valid in the linear regime of the dynamics. As the 
time increases the non-linear effect makes the situation different. We see structure even in 
the sample average distribution of density which develops slowly in a coherent manner. In 
Fig. 2 we plot {p{x,y;t)) at a time t = 150fm/c. This observed structure almost remains 
steady with increasing time, which corresponds to a equilibrium or at least a metastable 
state. The fluctuation is seen to be ~ 20% with respect to the initial uniform state. The 
appearance of such spatial inhomogeneity in the average density is in conformation with 



that of a mean trajectory calculation reported in Ref. |T2|. Though, the fluctuation in latter 
case is larger than what we observe here. This may be attributed to the different initial and 
also boundary conditions used in these two cases. However, large scale fluctuation exists in 
individual samples, as evident from the attainment of a non-zero steady value of in the 
present formalism. 

In conclusion, using a recently proposed simulation scheme for inclusion of fluctuation 
in BUU dynamics we study the problem of spinodal decomposition in nuclear matter. Our 
work is mainly devoted to the study of the large time behaviour of system. To analyize 
the simulation result we apply the well known Suzuki's model of dynamical phase transition 
successfully in the context of spinodal decomposition. Accordingly, the evolution of the 
system is guided by the co-operative effect of non-linearity and fluctuation in a universal 
manner. To make our stand clear, we would like to mention a recent work done by Baldo et 



al. [13 1 where the evolution of fluctuation is studied in the framework of Vlasov equation. In 
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the simulation work they studied the growth of an initial disturbance imposed on the uniform 
nuclear matter. Without being saturated, the overshoots very fast from its linear trends 
(in log scale). It will be relevant here to note that the finite collision rate, no matter how 
small it be, changes the nature of the solutions describing a diffusive process as modeled 
by the Langevin equation. Therefore, the non-diffusive behaviour of fluctuation, although 
exhibits some non-standard evolution pattern, but does not indicate conclusively that the 
dynamics of fragmentation is dominantly non-linear. However, similar such studies facilitate 



us to extract the growth constant tk, (the RPA frequency) of the initial state [0 and also 
allow a further scope to check the consistency of our work. The values of indicated in 



Fig. 1. provide a good fit to the dispersion relation shown in Fig 3. of Ref. ||10|| . 
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Figure Captions 

Fig.l The fluctuation coefficient L'^a^ versus time t for different values of unstable mode k 
{kx, ky) are shown. The open diamonds correspond to simulation results in steps of 
10 fm/c. The dashed curves represent the solution of Eq.(7) while solid curves show 
the nonlinear generalization of the same. The open circles represent the simulation 
result of a typical stable mode with [kx — 0.8 fm~^, ky = 0) which exhibits an early 
saturation. 

Fig. 2 The sample average density profile < p{x,y;t) > at time t= 150 fm/c versus the 
position {x, y) is shown. 
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